{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Richards equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#Simple 1D Richards' Equation Model\n", "\n", "This is a 1D Richards' equation model, written by Andrew Ireson, 3 November 2015. \n", "Downloaded from: https://github.com/amireson/RichardsEquation/tree/master \n", "Modified by Guo-Shiuan Lin, 2 August, 2023." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## In this assignment, you should get three different figures (1,2,3 as shown below) for each conditions (1,2,3,4). \n", "- If there are same plots (e.g., same soil properties from the same soil), don't need to show them in the report again.\n", "- Mainly show the plots where a change is expected to be observed" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [], "source": [ "# Import all of the basic libraries (you will always need these)\n", "from matplotlib import pyplot as pl\n", "import numpy as np\n", "\n", "# Import a library that contains soil moisture properties and functions\n", "# to import successfully, the vanGenuchten.py file must be in the same folder as this Jupyter notebook!\n", "import vanGenuchten as vg\n", "\n", "# Import ODE solvers\n", "from scipy.interpolate import interp1d\n", "from scipy.integrate import odeint\n", "\n", "# Select which soil properties to use. In the first example, we use the Hygiene Sand Stone. In the second example, you can try any soil you like\n", "p = vg.HygieneSandstone() " ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [], "source": [ "# Richards equation solver\n", "# This is a function that calculated the right hand side of Richards' equation. \n", "# You will not need to modify this function!\n", "# This block of code must be executed so that the function can be later called.\n", "\n", "def RichardsModel(psi,t,dz,n,p,vg,qTop,qBot,psiTop,psiBot):\n", " \n", " # Basic properties:\n", " C=vg.CFun(psi,p)\n", " \n", " # initialize vectors:\n", " q=np.zeros(n+1)\n", " \n", " # Upper boundary\n", " if qTop == []:\n", " KTop=vg.KFun(np.zeros(1)+psiTop,p)\n", " q[n]=-KTop*((psiTop-psi[n-1])/dz*2+1)\n", " else:\n", " q[n]=qTop\n", " \n", " # Lower boundary\n", " if qBot == []:\n", " if psiBot == []:\n", " # Free drainage\n", " KBot=vg.KFun(np.zeros(1)+psi[0],p)\n", " q[0]=-KBot\n", " else:\n", " # Type 1 boundary\n", " KBot=vg.KFun(np.zeros(1)+psiBot,p)\n", " q[0]=-KBot*((psi[0]-psiBot)/dz*2+1.0) \n", " else:\n", " # Type 2 boundary\n", " q[0]=qBot\n", " \n", " # Internal nodes\n", " i=np.arange(0,n-1)\n", " Knodes=vg.KFun(psi,p)\n", " Kmid=(Knodes[i+1]+Knodes[i])/2.0\n", " \n", " j=np.arange(1,n)\n", " q[j]=-Kmid*((psi[i+1]-psi[i])/dz+1.0)\n", " \n", " \n", " # Continuity\n", " i=np.arange(0,n)\n", " dpsidt=(-(q[i+1]-q[i])/dz)/C\n", " \n", " return dpsidt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Condition 1: use the default soil type" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Figure 1: hydraulic properties of the soil" ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(
,\n", " array([, , ], dtype=object))" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvwAAAEYCAYAAAAgZSKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAdZ0lEQVR4nO3cfZCV1X0H8N+yy+4mtrsW0BUDrpgmgHGaylJ5K2USdS1aZ/ijA06ngBmd6U6nVaS2gTATXyaTjW2TTpMI1mTRyYwaahDrTKlx/1AkQl9Cl0wbaE2jcUnKhllS76Kti8DpH3RvudkFuZd9PX4+M/ePezjPvefw7Hfnu8/efapSSikAAIAsTRrrBQAAACNH4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIyVXfhffvnluPXWW+Pyyy+PqqqqePbZZ9/zmF27dkVLS0vU19fHVVddFY888kglawVGiFxDfuQaGFB24X/77bfj4x//eHz1q189r/mvv/563HzzzbF06dLo6uqKz3zmM3HXXXfF9u3by14sMDLkGvIj18CAqpRSqvjgqqrYsWNHrFix4qxzPv3pT8dzzz0XBw8eLI61tbXF9773vdi7d2+lbw2MELmG/Mg1vL/VjPQb7N27N1pbW0vGbrrppujo6Ih33303Jk+ePOiY/v7+6O/vLz4/depU/OxnP4upU6dGVVXVSC8ZJryUUhw7diwuv/zymDRp+P9UR65h9Mk15Gekcz1gxAt/T09PNDU1lYw1NTXFiRMnore3N6ZPnz7omPb29njggQdGemmQvUOHDsWMGTOG/XXlGsaOXEN+RirXA0a88EfEoJ/yBz5FdLaf/jdu3Bjr168vPi8UCnHFFVfEoUOHoqGhYeQWCpno6+uLmTNnxi/+4i+O2HvINYwuuYb8jEauI0ah8F922WXR09NTMnbkyJGoqamJqVOnDnlMXV1d1NXVDRpvaGjwDQTKMFK/UpdrGDtyDfkZ6Y/Ajfh9+BctWhSdnZ0lYy+88ELMnz9/yM8DAuOfXEN+5BryVXbhf+utt2L//v2xf//+iDh9G6/9+/dHd3d3RJz+9d6aNWuK89va2uKNN96I9evXx8GDB2Pr1q3R0dER99577/DsALhgcg35kWugKJXpxRdfTBEx6LF27dqUUkpr165Ny5YtKznmpZdeStdee22qra1NV155ZdqyZUtZ71koFFJEpEKhUO5y4X2p3MzINYx/cg35Ga3MXNB9+EdLX19fNDY2RqFQ8JlAOA8TITMTYY0wnkyEzEyENcJ4MlqZGfHP8AMAAGNH4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIwp/AAAkLGKCv/mzZtj1qxZUV9fHy0tLbF79+5zzn/iiSfi4x//eHzwgx+M6dOnx6c+9ak4evRoRQsGRoZcQ37kGoiooPBv27Yt1q1bF5s2bYqurq5YunRpLF++PLq7u4ec/53vfCfWrFkTd9xxR3z/+9+Pp59+Ov7pn/4p7rzzzgtePDA85BryI9dAUSrTddddl9ra2krG5syZkzZs2DDk/D/7sz9LV111VcnYl7/85TRjxozzfs9CoZAiIhUKhXKXC+9L5WZGrmH8k2vIz2hlpqwr/MePH499+/ZFa2tryXhra2vs2bNnyGMWL14cP/7xj2Pnzp2RUoqf/vSn8a1vfStuueWWs75Pf39/9PX1lTyAkSHXkB+5Bs5UVuHv7e2NkydPRlNTU8l4U1NT9PT0DHnM4sWL44knnohVq1ZFbW1tXHbZZXHxxRfHV77ylbO+T3t7ezQ2NhYfM2fOLGeZQBnkGvIj18CZKvqj3aqqqpLnKaVBYwMOHDgQd911V3z2s5+Nffv2xfPPPx+vv/56tLW1nfX1N27cGIVCofg4dOhQJcsEyiDXkB+5BiIiasqZPG3atKiurh50deDIkSODriIMaG9vjyVLlsQf//EfR0TEr/zKr8RFF10US5cujc997nMxffr0QcfU1dVFXV1dOUsDKiTXkB+5Bs5U1hX+2traaGlpic7OzpLxzs7OWLx48ZDH/Pd//3dMmlT6NtXV1RFx+koDMLbkGvIj18CZyv5Iz/r16+PrX/96bN26NQ4ePBj33HNPdHd3F3/lt3HjxlizZk1x/q233hrPPPNMbNmyJV577bV45ZVX4q677orrrrsuLr/88uHbCVAxuYb8yDUwoKyP9ERErFq1Ko4ePRoPPvhgHD58OK655prYuXNnNDc3R0TE4cOHS+7xe/vtt8exY8fiq1/9avzRH/1RXHzxxfHJT34yHnrooeHbBXBB5BryI9fAgKo0AX5P19fXF42NjVEoFKKhoWGslwPj3kTIzERYI4wnEyEzE2GNMJ6MVmYquksPAAAwMSj8AACQMYUfAAAypvADAEDGFH4AAMiYwg8AABlT+AEAIGMKPwAAZEzhBwCAjCn8AACQMYUfAAAypvADAEDGFH4AAMiYwg8AABlT+AEAIGMKPwAAZEzhBwCAjCn8AACQMYUfAAAypvADAEDGFH4AAMiYwg8AABlT+AEAIGMKPwAAZEzhBwCAjCn8AACQMYUfAAAypvADAEDGFH4AAMiYwg8AABlT+AEAIGMKPwAAZEzhBwCAjCn8AACQMYUfAAAyVlHh37x5c8yaNSvq6+ujpaUldu/efc75/f39sWnTpmhubo66urr48Ic/HFu3bq1owcDIkGvIj1wDERE15R6wbdu2WLduXWzevDmWLFkSf/VXfxXLly+PAwcOxBVXXDHkMStXroyf/vSn0dHREb/8y78cR44ciRMnTlzw4oHhIdeQH7kGBlSllFI5ByxYsCDmzZsXW7ZsKY7NnTs3VqxYEe3t7YPmP//883HbbbfFa6+9FlOmTKlokX19fdHY2BiFQiEaGhoqeg14Pyk3M3IN459cQ35GKzNlfaTn+PHjsW/fvmhtbS0Zb21tjT179gx5zHPPPRfz58+PP/3TP40PfehD8dGPfjTuvffe+J//+Z/KVw0MG7mG/Mg1cKayPtLT29sbJ0+ejKamppLxpqam6OnpGfKY1157Lb7zne9EfX197NixI3p7e+P3f//342c/+9lZPxfY398f/f39xed9fX3lLBMog1xDfuQaOFNFf7RbVVVV8jylNGhswKlTp6KqqiqeeOKJuO666+Lmm2+OL33pS/H444+f9apBe3t7NDY2Fh8zZ86sZJlAGeQa8iPXQESZhX/atGlRXV096OrAkSNHBl1FGDB9+vT40Ic+FI2NjcWxuXPnRkopfvzjHw95zMaNG6NQKBQfhw4dKmeZQBnkGvIj18CZyir8tbW10dLSEp2dnSXjnZ2dsXjx4iGPWbJkSfznf/5nvPXWW8WxV199NSZNmhQzZswY8pi6urpoaGgoeQAjQ64hP3INlEhl+uY3v5kmT56cOjo60oEDB9K6devSRRddlH70ox+llFLasGFDWr16dXH+sWPH0owZM9Jv//Zvp+9///tp165d6SMf+Ui68847z/s9C4VCiohUKBTKXS68L5WbGbmG8U+uIT+jlZmy78O/atWqOHr0aDz44INx+PDhuOaaa2Lnzp3R3NwcERGHDx+O7u7u4vxf+IVfiM7OzvjDP/zDmD9/fkydOjVWrlwZn/vc54bj5xVgGMg15EeugQFl34d/LLivL5RnImRmIqwRxpOJkJmJsEYYT8blffgBAICJReEHAICMKfwAAJAxhR8AADKm8AMAQMYUfgAAyJjCDwAAGVP4AQAgYwo/AABkTOEHAICMKfwAAJAxhR8AADKm8AMAQMYUfgAAyJjCDwAAGVP4AQAgYwo/AABkTOEHAICMKfwAAJAxhR8AADKm8AMAQMYUfgAAyJjCDwAAGVP4AQAgYwo/AABkTOEHAICMKfwAAJAxhR8AADKm8AMAQMYUfgAAyJjCDwAAGVP4AQAgYwo/AABkTOEHAICMKfwAAJAxhR8AADJWUeHfvHlzzJo1K+rr66OlpSV27959Xse98sorUVNTE7/6q79aydsCI0iuIT9yDURUUPi3bdsW69ati02bNkVXV1csXbo0li9fHt3d3ec8rlAoxJo1a+L666+veLHAyJBryI9cAwOqUkqpnAMWLFgQ8+bNiy1bthTH5s6dGytWrIj29vazHnfbbbfFRz7ykaiuro5nn3029u/ff97v2dfXF42NjVEoFKKhoaGc5cL7UrmZkWsY/+Qa8jNamSnrCv/x48dj37590draWjLe2toae/bsOetxjz32WPzwhz+M++6777zep7+/P/r6+koewMiQa8iPXANnKqvw9/b2xsmTJ6OpqalkvKmpKXp6eoY85gc/+EFs2LAhnnjiiaipqTmv92lvb4/GxsbiY+bMmeUsEyiDXEN+5Bo4U0V/tFtVVVXyPKU0aCwi4uTJk/E7v/M78cADD8RHP/rR8379jRs3RqFQKD4OHTpUyTKBMsg15EeugYiI8/sR/v9MmzYtqqurB10dOHLkyKCrCBERx44di+9+97vR1dUVf/AHfxAREadOnYqUUtTU1MQLL7wQn/zkJwcdV1dXF3V1deUsDaiQXEN+5Bo4U1lX+Gtra6OlpSU6OztLxjs7O2Px4sWD5jc0NMS//Mu/xP79+4uPtra2mD17duzfvz8WLFhwYasHLphcQ37kGjhTWVf4IyLWr18fq1evjvnz58eiRYvi0Ucfje7u7mhra4uI07/e+8lPfhLf+MY3YtKkSXHNNdeUHH/ppZdGfX39oHFg7Mg15EeugQFlF/5Vq1bF0aNH48EHH4zDhw/HNddcEzt37ozm5uaIiDh8+PB73uMXGF/kGvIj18CAsu/DPxbc1xfKMxEyMxHWCOPJRMjMRFgjjCfj8j78AADAxKLwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADImMIPAAAZU/gBACBjCj8AAGRM4QcAgIwp/AAAkDGFHwAAMqbwAwBAxhR+AADIWEWFf/PmzTFr1qyor6+PlpaW2L1791nnPvPMM3HjjTfGJZdcEg0NDbFo0aL49re/XfGCgZEh15AfuQYiKij827Zti3Xr1sWmTZuiq6srli5dGsuXL4/u7u4h57/88stx4403xs6dO2Pfvn3xiU98Im699dbo6uq64MUDw0OuIT9yDQyoSimlcg5YsGBBzJs3L7Zs2VIcmzt3bqxYsSLa29vP6zU+9rGPxapVq+Kzn/3sec3v6+uLxsbGKBQK0dDQUM5y4X2p3MzINYx/cg35Ga3MlHWF//jx47Fv375obW0tGW9tbY09e/ac12ucOnUqjh07FlOmTDnrnP7+/ujr6yt5ACNDriE/cg2cqazC39vbGydPnoympqaS8aampujp6Tmv1/jiF78Yb7/9dqxcufKsc9rb26OxsbH4mDlzZjnLBMog15AfuQbOVNEf7VZVVZU8TykNGhvKU089Fffff39s27YtLr300rPO27hxYxQKheLj0KFDlSwTKINcQ37kGoiIqCln8rRp06K6unrQ1YEjR44Muorw87Zt2xZ33HFHPP3003HDDTecc25dXV3U1dWVszSgQnIN+ZFr4ExlXeGvra2NlpaW6OzsLBnv7OyMxYsXn/W4p556Km6//fZ48skn45ZbbqlspcCIkGvIj1wDZyrrCn9ExPr162P16tUxf/78WLRoUTz66KPR3d0dbW1tEXH613s/+clP4hvf+EZEnP7msWbNmvjLv/zLWLhwYfFqwwc+8IFobGwcxq0AlZJryI9cA0WpAg8//HBqbm5OtbW1ad68eWnXrl3Ff1u7dm1atmxZ8fmyZctSRAx6rF279rzfr1AopIhIhUKhkuXC+04lmZFrGN/kGvIzWpkp+z78Y8F9faE8EyEzE2GNMJ5MhMxMhDXCeDIu78MPAABMLAo/AABkTOEHAICMKfwAAJAxhR8AADKm8AMAQMYUfgAAyJjCDwAAGVP4AQAgYwo/AABkTOEHAICMKfwAAJAxhR8AADKm8AMAQMYUfgAAyJjCDwAAGVP4AQAgYwo/AABkTOEHAICMKfwAAJAxhR8AADKm8AMAQMYUfgAAyJjCDwAAGVP4AQAgYwo/AABkTOEHAICMKfwAAJAxhR8AADKm8AMAQMYUfgAAyJjCDwAAGVP4AQAgYwo/AABkTOEHAICMVVT4N2/eHLNmzYr6+vpoaWmJ3bt3n3P+rl27oqWlJerr6+Oqq66KRx55pKLFAiNHriE/cg1EVFD4t23bFuvWrYtNmzZFV1dXLF26NJYvXx7d3d1Dzn/99dfj5ptvjqVLl0ZXV1d85jOfibvuuiu2b99+wYsHhodcQ37kGhhQlVJK5RywYMGCmDdvXmzZsqU4Nnfu3FixYkW0t7cPmv/pT386nnvuuTh48GBxrK2tLb73ve/F3r17z+s9+/r6orGxMQqFQjQ0NJSzXHhfKjczcg3jn1xDfkYrMzXlTD5+/Hjs27cvNmzYUDLe2toae/bsGfKYvXv3Rmtra8nYTTfdFB0dHfHuu+/G5MmTBx3T398f/f39xeeFQiEiTv+nAO9tICvn8/O8XMPEINeQn3JyfSHKKvy9vb1x8uTJaGpqKhlvamqKnp6eIY/p6ekZcv6JEyeit7c3pk+fPuiY9vb2eOCBBwaNz5w5s5zlwvve0aNHo7Gx8Zxz5BomFrmG/JxPri9EWYV/QFVVVcnzlNKgsfeaP9T4gI0bN8b69euLz998881obm6O7u7uEf3PGA19fX0xc+bMOHTo0IT/dWdOe4nIaz+FQiGuuOKKmDJlynkfI9eVy+lrJyKv/eS0F7keXTl97UTktZ+c9lJJritRVuGfNm1aVFdXD7o6cOTIkUFXBQZcdtllQ86vqamJqVOnDnlMXV1d1NXVDRpvbGyc8Cd2QENDg72MUzntZ9Kk9/67fLkePjl97UTktZ+c9iLXoyunr52IvPaT017OJ9cX9PrlTK6trY2Wlpbo7OwsGe/s7IzFixcPecyiRYsGzX/hhRdi/vz5Q34eEBhdcg35kWvgTGX/OLF+/fr4+te/Hlu3bo2DBw/GPffcE93d3dHW1hYRp3+9t2bNmuL8tra2eOONN2L9+vVx8ODB2Lp1a3R0dMS99947fLsALohcQ37kGihKFXj44YdTc3Nzqq2tTfPmzUu7du0q/tvatWvTsmXLSua/9NJL6dprr021tbXpyiuvTFu2bCnr/d5555103333pXfeeaeS5Y4r9jJ+5bSfSvYi15XLaS8p5bWf9/te5LpyOe0lpbz2Yy/lK/s+/AAAwMQxsn8hAAAAjCmFHwAAMqbwAwBAxhR+AADI2JgU/s2bN8esWbOivr4+WlpaYvfu3eecv2vXrmhpaYn6+vq46qqr4pFHHhk0Z/v27XH11VdHXV1dXH311bFjx46RWv4g5eznmWeeiRtvvDEuueSSaGhoiEWLFsW3v/3tkjmPP/54VFVVDXq88847I72Vsvby0ksvDbnOf/u3fyuZN1bnppy93H777UPu5WMf+1hxzlidl5dffjluvfXWuPzyy6OqqiqeffbZ9zxmLDIj13I9GuR69DOTU7blWq5H0rjO9YjeA2gI3/zmN9PkyZPT1772tXTgwIF09913p4suuii98cYbQ85/7bXX0gc/+MF09913pwMHDqSvfe1rafLkyelb3/pWcc6ePXtSdXV1+vznP58OHjyYPv/5z6eampr093//9+NuP3fffXd66KGH0j/+4z+mV199NW3cuDFNnjw5/fM//3NxzmOPPZYaGhrS4cOHSx7jbS8vvvhiioj07//+7yXrPHHiRHHOWJ2bcvfy5ptvluzh0KFDacqUKem+++4rzhmr87Jz5860adOmtH379hQRaceOHeecPxaZkWu5luvyTIRcp5RXtuVarkfaeM71qBf+6667LrW1tZWMzZkzJ23YsGHI+X/yJ3+S5syZUzL2e7/3e2nhwoXF5ytXrky/+Zu/WTLnpptuSrfddtswrfrsyt3PUK6++ur0wAMPFJ8/9thjqbGxcbiWeN7K3cvAN5D/+q//OutrjtW5udDzsmPHjlRVVZV+9KMfFcfG6ryc6Xy+gYxFZuR6MLkefnI9+pnJKdtyLdejabzlelQ/0nP8+PHYt29ftLa2loy3trbGnj17hjxm7969g+bfdNNN8d3vfjfefffdc84522sOl0r28/NOnToVx44diylTppSMv/XWW9Hc3BwzZsyI3/qt34qurq5hW/dQLmQv1157bUyfPj2uv/76ePHFF0v+bSzOzXCcl46Ojrjhhhuiubm5ZHy0z0slRjszcj2YXA8/uR79zOSUbbk+Ta7Hl9HMy6gW/t7e3jh58mQ0NTWVjDc1NUVPT8+Qx/T09Aw5/8SJE9Hb23vOOWd7zeFSyX5+3he/+MV4++23Y+XKlcWxOXPmxOOPPx7PPfdcPPXUU1FfXx9LliyJH/zgB8O6/jNVspfp06fHo48+Gtu3b49nnnkmZs+eHddff328/PLLxTljcW4u9LwcPnw4/u7v/i7uvPPOkvGxOC+VGO3MyPVgcj385Hr0M5NTtuVart/vua65sKVWpqqqquR5SmnQ2HvN//nxcl9zOFX63k899VTcf//98Td/8zdx6aWXFscXLlwYCxcuLD5fsmRJzJs3L77yla/El7/85eFb+BDK2cvs2bNj9uzZxeeLFi2KQ4cOxZ//+Z/Hb/zGb1T0msOp0vd9/PHH4+KLL44VK1aUjI/leSnXWGRGrk+T65El1/9vtDKTU7bl+jS5Hj9GKy+jeoV/2rRpUV1dPeinkiNHjgz66WXAZZddNuT8mpqamDp16jnnnO01h0sl+xmwbdu2uOOOO+Kv//qv44Ybbjjn3EmTJsWv/dqvjehPpheylzMtXLiwZJ1jcW4uZC8ppdi6dWusXr06amtrzzl3NM5LJUY7M3L9/+RarkfKWGQmp2zL9WByPfZGMy+jWvhra2ujpaUlOjs7S8Y7Oztj8eLFQx6zaNGiQfNfeOGFmD9/fkyePPmcc872msOlkv1EnL5ScPvtt8eTTz4Zt9xyy3u+T0op9u/fH9OnT7/gNZ9NpXv5eV1dXSXrHItzcyF72bVrV/zHf/xH3HHHHe/5PqNxXiox2pmR69PkWq5H0lhkJqdsy/Vgcj32RjUvZf2J7zAYuP1SR0dHOnDgQFq3bl266KKLin9dvWHDhrR69eri/IFbFt1zzz3pwIEDqaOjY9Ati1555ZVUXV2dvvCFL6SDBw+mL3zhC6N++77z3c+TTz6Zampq0sMPP1xyq6g333yzOOf+++9Pzz//fPrhD3+Yurq60qc+9alUU1OT/uEf/mFc7eUv/uIv0o4dO9Krr76a/vVf/zVt2LAhRUTavn17cc5YnZty9zLgd3/3d9OCBQuGfM2xOi/Hjh1LXV1dqaurK0VE+tKXvpS6urqKtywbD5mRa7mW6/JMhFynlFe25Vqu38+5HvXCn1JKDz/8cGpubk61tbVp3rx5adeuXcV/W7t2bVq2bFnJ/Jdeeilde+21qba2Nl155ZVpy5Ytg17z6aefTrNnz06TJ09Oc+bMKfkiHmnl7GfZsmUpIgY91q5dW5yzbt26dMUVV6Ta2tp0ySWXpNbW1rRnz55xt5eHHnooffjDH0719fXpl37pl9Kv//qvp7/9278d9JpjdW7K/Tp788030wc+8IH06KOPDvl6Y3VeBm6ndravmfGSGbmW69Eg16OfmZyyLddyPZLGc66rUvq/vw4AAACyM6qf4QcAAEaXwg8AABlT+AEAIGMKPwAAZEzhBwCAjCn8AACQMYUfAAAypvADAEDGFH4AAMiYwg8AABlT+AEAIGMKPwAAZOx/ATPSJvdPqgpSAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Three plots. Plot three soil properties (theta, C, K calculated below) on y with psi on x. \n", "psi = np.linspace(-10,0)\n", "theta = vg.thetaFun(psi,p)\n", "C=vg.CFun(psi,p)\n", "K=vg.KFun(psi,p)\n", "# Your code to continue here. Use pl.subplots(1,3, sharex=True) to create 3 plots \n", "pl.subplots(1,3, sharex=True, figsize=(9,3))" ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\glin\\AppData\\Local\\Temp\\ipykernel_39692\\4283632028.py:30: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)\n", " q[0]=-KBot*((psi[0]-psiBot)/dz*2+1.0)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Model run successfully\n" ] } ], "source": [ "# This block of code sets up and runs the model. You must run it again after you change some parameters \n", "\n", "# Boundary conditions\n", "qTop=-0.01 # m/day. It is negative because the direction of infiltration is downward\n", "qBot=[]\n", "psiTop=[] # top boundary condition \n", "psiBot= 0 # set to [] to change to set bottom boundary condition as free drainage \n", "\n", "# Grid in space\n", "dz=0.1 # meter\n", "ProfileDepth=5 # meter\n", "z=np.arange(dz/2.0,ProfileDepth,dz)\n", "n=z.size\n", "\n", "# Grid in time\n", "t = np.linspace(0,10,100) # 10 days with 100 time steps\n", "\n", "# Initial conditions\n", "psi0=-z\n", "\n", "# Solve. The odeint function solves the system\n", "psi=odeint(RichardsModel,psi0,t,args=(dz,n,p,vg,qTop,qBot,psiTop,psiBot),mxstep=5000000);\n", "\n", "print(\"Model run successfully\") " ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [], "source": [ "# Post process model output to get some variables of interest\n", "\n", "# Get water content\n", "theta=vg.thetaFun(psi,p)\n", "\n", "# Get total profile storage\n", "S=theta.sum(axis=1)*dz\n", "\n", "# Get change in storage [dVol]\n", "dS=np.zeros(S.size)\n", "dS[1:]=np.diff(S)/(t[1]-t[0])\n", "\n", "# Get infiltration flux\n", "if qTop == []:\n", " KTop=vg.KFun(np.zeros(1)+psiTop,p)\n", " qI=-KTop*((psiTop-psi[:,n-1])/dz*2+1)\n", "else:\n", " qI=np.zeros(t.size)+qTop\n", " \n", "# Get discharge flux\n", "if qBot == []:\n", " if psiBot == []:\n", " # Free drainage\n", " KBot=vg.KFun(psi[:,0],p)\n", " qD=-KBot\n", " else:\n", " # Type 1 boundary\n", " KBot=vg.KFun(np.zeros(1)+psiBot,p)\n", " qD=-KBot*((psi[:,0]-psiBot)/dz*2+1.0)\n", "else:\n", " qD=np.zeros(t.size)+qBot\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Figure 2: Vertical profile of the soil" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [], "source": [ "## Plot vertical profiles. Uncomment to see the plot\n", "\n", "# pl.rcParams['figure.figsize'] = (10.0, 10.0)\n", "# for i in range(0,t.size,1): ### set the time to plot\n", "# pl.subplot(121)\n", "# pl.plot(psi[i,:],z)\n", "# pl.subplot(122)\n", "# pl.plot(theta[i,:],z)\n", "\n", "# pl.subplot(121)\n", "# pl.ylabel('Height [m]',fontsize=20) \n", "# # the Height for soil may seem confusing. Imaging there is a soil column standing on the ground. \n", "# # The bottom of the soil colum touching the ground is height=0\n", "# # the top of the soil column 5 meters above the ground is height=5\n", "# pl.xlabel(r'$\\psi$ [m]',fontsize=20)\n", "# pl.subplot(122)\n", "# pl.xlabel(r'$\\theta$ [-]',fontsize=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Figure 3: temporal change of the water flows" ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [], "source": [ "## Plot timeseries, Uncomment to see the plot\n", "\n", "# dt = t[2]-t[1]\n", "# pl.plot(t,dS,label='Change in storage')\n", "# pl.plot(t,-qI,label='Infiltration')\n", "# pl.plot(t,-qD,label='Discharge')\n", "# pl.legend(loc=1)\n", "# pl.xlabel('days')\n", "# pl.ylabel('$m^3$/days')" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Condition 2. Change the soil type and compare the results\n", "- You may want to change the simuulation duration t to see the change in discharge" ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [], "source": [ "# Choose another soil type\n", "p=vg.GuelphLoamDrying() # or p=vg.BeitNetofaClay(), p=vg.GuelphLoamWetting(), p=vg.SiltLoamGE3(), p=vg.TouchetSiltLoam() \n", "\n", "# run the model and plot again" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Condition 3. With the same soil, change the infiltration rate and compare\n", "- You may want to change the simuulation duration t to shorter (if you increase qTop) or longer (if you reduce qTop) to see the change in discharge" ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [], "source": [ "qTop=-0.001 # or -0.005 m/day or other number you want to try\n", "\n", "# run the model and plot again" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Condition 4. With the same soil, change the boundary condition to free drainage and compare" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [], "source": [ "# run the model and plot again" ] } ], "metadata": { "kernelspec": { "display_name": "base", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 1 }